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ABSTRACT 

An automatic history matching algorithm is developed based on bi-cubic 

spline approximations of permeability and porosity in a single-phase, two- 

dimensional areal reservoir from well pressure data. The regularization 

feature of the algorithm, the theoretical details of which are described by 

Kravaris and Seinfeld^’ ^ is used to convert the ill-posed history matching 

problem into a well-posed problem. The algorithm employs the conjugate 

g 

gradient method of Nazareth as its core minimization method. A number of 
numerical experiments are carried out to evaluate the performance of the 
algorithm. Comparisons with conventional (non-regularized) automatic history 
matching algorithms Indicate the superiority of the new algorithm with respect 
to the parameter estimates obtained. A quasioptimal regularization parameter 
is determined without requiring a priori information on the statistical 
properties of the observations. 


The research for the third author was partially supported under the National 
Aeronautics and Space Administration under NASA Contract Nos. NAS1-17070 and 
NAS1-18107 while he was in residence at the Institute for Computer 
Applications in Science and Engineering (ICASE), NASA Langley Research Center, 
Hampton, VA 23665. 
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INTRODUCTION 


The process of estimating unknown properties, such as permeability and 
porosity, in a mathematical reservoir model to give the best fit to measured 
well pressure and production data is commonly called "history matching." 
Because the properties in an inhomogeneous reservoir vary with location, 
conceptually an infinite number of parameters is required for a full 
description of the reservoir. From a computational point of view, a reservoir 
simulator contains only a finite number of parameters, corresponding to each 
grid block or element in the spatial domain. In field scale simulations, it is 
not unusual for the reservoir domain to consist of the order of 10,000 grid 
blocks, and consequently 20,000 or more parameters may need to be estimated 
simultaneously. This potential large dimensionality of the unknown parameters 
distinguishes the reservoir history matching problem from other parameter 
estimation problems in science and engineering. Moreover, the standard 
reservoir history matching problem is mathematically ill- posed, and this 
inherent ill-posedness, coupled with such a large number of unknown 
parameters, lies at the root of the difficulties in its solution. The ill- 
posedness of the history matching problem manifests itself by numerical 
instabilities in the estimated parameters. Such instabilities are well 
documented in the petroleum engineering and hydrology literature 1 * 2 . 

The principal approach that has been used to alleviate ill- conditioning in 
the parameter estimates is to decrease the number of unknown parameters, and, 
if possible, utilize any available information to constrain the space of 
unknown parameters. One commonly used approach for reducing the number of 
unknown parameters is to divide the reservoir into a relatively small number 
of zones and to assume uniform properties within each zone. While this 
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approach is effective in reducing the number of unknowns, sufficient a priori 
information is usually not available to enable specification of the zones on 
any physical basis. Moreover, even though limiting the dimension of the 
parameter space, zonation does not alleviate the fundamental ill- posed nature 
of the problem. An alternative to zonation is to use prior information 
expressed in the form of an assumed probability distribution for the reservoir 
properties. If certain a priori knowledge is assumed about the parameter mean 
values and correlations, the history matching objective function can be 
modified to include a term penalizing the weighted deviations of the 
parameters from their assumed mean values 1 . A form of Bayesian estimation can 
then be used to determine the unknown parameters. While it has been shown 
that better conditioned estimates may be obtained when a priori statistical 
information is employed than when it is not sufficient knowledge of the nature 
of the unknown parameters is generally not available to specify the parameters 
needed to carry out a Bayesian estimation. 

The critical problems in generating an effective algorithm for history 
matching are twofold: (1) The original problem must be defined in a manner 
that alleviates the ill-posed nature of the problem; and (2) An efficient 
computational algorithm must be developed for solving the large, constrained, 
nonlinear minimization problem that results from any history matching problem. 

With respect to the inherent ill-posedness of the history matching 
problem, Kravaris and Seinfeld 3 * 4 have shown that the concept of regularization 
can be applied rigorously to the estimation of spatially-varying parameters in 
partial differential equations of parabolic type. The regularization ’idea, 
first advanced by Tikhonov 18 , has been widely used in the solution of ill-posed 
integral equations, but had not heretofore been developed for the estimation 



3 


of parameters in partial differential equations. In short, regularization of a 
problem refers to solving a related problem, called the regularized problem, 
whose solution is more "regular" (in a certain sense) than the solution of the 
original problem and approximates (in a certain sense) the solution of the 
original problem. More precisely, regularization of an ill- posed problem 
refers to solving a well- posed problem, whose solution gives a physically 
meaningful answer to the original ill-posed problem. The regularization 
formulation of parameter estimation measures the "non-smoothness" of the 
estimated parameter as a norm of the parameter in an appropriate Hilbert 
space. No prior information about the parameter is required other than a 
general idea of the degree of smoothness desired in the estimated field. The 
only unspecified parameter is that reflecting the relative weight given to the 
smoothness norm versus the usual least-square objective function. 

Banks and co-workers 5-8 have found that the use of spline representations 
for spatially- varying parameters in one-dimensional partial differential 
equations of both parabolic and hyperbolic type leads to well-conditioned 
estimates. Although their numerical results were obtained for low levels of 
spline discretization, it seems that the spline representation itself may 
impart a degree of smoothness to the parameter distribution that could 
circumvent some of the ill- conditioning inherent in the finite-difference or 
zonation representation of the unknown parameters. The use of two- 
dimensional, bi-cubic spline approximations for reservoir history matching is 
an additional new feature of the work reported here. 

The object of this work is to present an automatic history matching 
algorithm based on the concept of regularization together with bi-cubic spline 
approximations for the estimation of permeability and porosity in a single- 
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phase, two-dimensional areal reservoir simulation. The two critical problems 
cited above are addressed in the algorithm. First, the regularization 
formulation converts the history matching problem to a mathematically well- 
posed problem. Second, we use a particularly efficient numerical minimization 
method, the conjugate gradient method of Nazareth 9 , as the core technique for 
the minimization. We present the results of extensive numerical testing of the 
algorithm in which both permeability and porosity distributions are estimated. 
Of particular interest will be the effects of the degree of regularization and 
of the order of the spline approximation on the behavior of the estimates. 
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HISTORY MATCHING BY REGULARIZATION 

The problem of history matching may be viewed in a general way by 
expressing the reservoir model, or simulator, as the nonlinear operator 
equation, 

Ka = u (1) 

where a represents the reservoir parameters, K is the operator representing 
the reservoir model, and u is the observed portion of the model’s output, such 
as the well pressures. The history matching problem is just the inverse 
problem to Eq. (1), that is given u and K find a. This inverse problem is 
well-posed if: 

(a) For every u there exists a solution a 

(b) The solution is unique 

(c) The solution is stable, i.e., small perturbations in u imply small 
perturbations in a. 

If any of (a), (b) or (c) is false, the inverse problem is ill-posed. 
Establishing uniqueness of a given u for operators K typical of reservoir 
simulators is an extremely difficult problem, and at this time uniqueness 
results are only available for very special cases 10 . It can be shown readily, 
however, that the inverse problem to Eq. (1) for spatially- varying parameters 
in parabolic partial differential equations is unstable 3 , and the estimation of 
a from u is an ill-posed problem. As we noted in the Introduction, the ill- 
posedness manifests itself by highly ill-conditioned estimates in conventional 
automatic history matching approaches. 

Let us now be more specific and consider unsteady flow of a slightly 
compressible oil with viscosity p in a two dimensional, areal reservoir of unit 
thickness, spatial domain ft and boundary 9ft in which fluid is being withdrawn 
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from N w wells located at (x w ,y w ), w = 1,2 N w . The fluid properties, p and 

c, are assumed to be known whereas the porosity <f> and permeability k are 
assumed to be unknown. The pressure distribution in the reservoir is governed 
by 

N w 

c<fl It = v * ( u Vp ^ + 21 q w S(x-x w )6(y-y w ) (2) 

w=1 

in ft x[0,T] 

|£ - 0 on an x [0,T] (3) 

p(x,y,0) = p 0 (x,y) in n (4) 

where 8p/an is the outward normal derivative of p on the boundary an and [0,T] 
is the time interval over which data are available. Because of the small size 
of the well bores relative to the reservoir dimensions, the well flow rates 
are represented as point sink terms in the pressure equation. If there exist 
observed pressures at Nt times at N 0 bs locations, Pk,n ot)S » k ■ 1,2,..., N 0 t> s , n 
= 1,2,..., Nt, then the customary history matching least-squares objective 
function is 

N t N obs 

J LS = Z Z (Pk,n obs ~ P( x k>yk>t n )) 2 (5) 

n=1 k=1 

The conventional history matching problem can be viewed therefore as a non- 
linear optimization problem of minimizing the sum of squares of differences 
between the observed and predicted pressures subject to the constraint of the 
reservoir model, Eqs. (2)-(-4) . 
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In the regularization approach we minimize an augmented objective 
function, called the smoothing functional, denoted by Jsm> that consists of the 
sum of the least-squares term, Jls, and a stabilizing functional, Jgj. The 
stabilizing functional for a parameter a(a=<{> or k) is of the form 

J ST = I l a l | 2 H 3 (ft) 

where | | | | 2 h 3 (q) is a norm defined in the Sobolev space H 3 (ft).* Thus the 
overall objective function to be minimized is 

J SM = J LS + &a J ST (7) 

Where Ba is a weighting coefficient chosen to reflect the degree of importance 
given to Jg-p. 

The minimization of Jg^ is performed over an appropriate finite- 
dimensional subspace of H 3 (fi), the so-called space of approximants which can be 
spanned by cubic spline functions. Thus the infinite-dimensional parameter 
spaces for k and <(> are reduced to finite-dimensional spaces by cubic spline 
approximations, and the finite dimensional minimization of Jg^ is carried out 
by an appropriate numerical minimization method. 


*The Sobolev space H 3 (G) is the set of functions that are square- integrable 
over Q and have square- integrable derivatives up to order 3- The norm of H 3 (fi) 
is given by equation (A. 12) of Appendix A. 
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BICUBIC SPLINE APPROXIMATION OF PERMEABILITY AND POROSITY 
A general approach to representing the spatial variation of reservoir 
properties is through the use of bicubic spline functions, in which a parameter 
a(x,y) is represented as 


a(x,y) 


Nys 

z 



b x (ix> x )W a 5, Xt iy byUy.y) 


( 8 ) 


where 

b x^x> x ) = X*^~*x + *x = 1.2,..., N xs (9) 

by(%,y) = X^^-fcy + a y = 1,2 Nys (10) 


and where x***( ) is the cubic B-spline function, 
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x 3 

T 


x e [0,1] 


1 . x-1 . (x-1) 2 (x— 1 ) 3 

F ~ 2 ~ 2 ' 


x e [1,2] 


x* 1, Cx) - 



+ 


(x-2) 3 

2 


x e [2,3] 


1 _ x-3 . (x-3) 2 _ (x-3) 3 

F ~ ~ 2 ~ ~ 2 ~ ~ — F~ 


x e [3, 4] 


0 


otherwise 


( 11 ) 


Ax s and Ay s are the grid distances of the spline grid along the x- and in- 
directions, respectively. With this approximation, a(x,y) is replaced by the 
set of unknown coefficients, W a £ Xf £ , X, x = 1,2,..., N xs and £y = 1,2,..., Ny S . 

The grid used for spline representation of the unknown properties need 
not necessarily coincide with that on which the actual reservoir model is 
solved. The reservoir model will be solved numerically using finite-difference 
approximations on the block-centered grid system shown in Figure 1. Figure 1 
also shows the spline grid system. The finite-difference grid can be expressed 
compactly as N x = {i x |i x = 1,2,..., N x }, Ny = { iy | iy = 1,2,..., Ny} and % = {i|i 
= i^ + N x (iy— 1), i x e N x , iy e Ny} = = 1 »2, •••, N}, where N = N x Ny. 

The finite-difference approximation of the pressure equation can be 
written in compact notation as 
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- W n->) - Z QL,j ^ (Pj n - Pi") * I q w «i,i„ 

jeJi w=1 


where 


for i e Njj, n = 1,2 Nt 


Ji - (j|j = i“N x , i-1 , i+1 , i+N x } D Njj 
Q = AxAy/At 


Ay/ Ax if j = i-1, i+1 and j e % 


Ax/ Ay if j = i~N x , i+N x and j e % 


6 i,i w is the Kroneeker delta, and k(i,j) = (ki + kj)/2. The initial condition is 
Pi 0 = Po» i e %. 

The least-squares objective function is then written as 


N N t N obs 


•JLS - Z r Z (Pk,n obs - Pi n > «i,i k 

* ■ 4 A I . -t •» 


1=1 n=1 k=1 


where we assume that t n - t n _-| = At for n = 1, ..., Nt* 
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HISTORY MATCHING ALGORITHM 

The problem we now seek to solve is to minimize the augmented objective 
function Jg^, given by Eq. (7), with respect to the spline coefficients 

W a £ x o , £ x = 1,2 N xs and £y = 1,2,..., Ny S , subject to the pressure 

equation (12). To obtain an algorithm to solve this problem two steps are 
required. First, we must compute the gradient of Jg^ with respect to each 
W a £ x ,j,yf and, second, that gradient is then used in a numerical minimization 
method to minimize JgM* The calculation of these gradients represents the 
most time consuming part of updating the parameter iterates. In a problem as 
large as history matching in a candidate algorithm these derivatives must be 
able to be calculated directly, not requiring the individual derivatives 
3Pi n /3W a o . . 11 » 15 Those algorithms based on an optimal control, or 

x ,Jty 

variational, formulation possess this necessary property. First, we solve the 
reservoir simulator equation, Eq. (13), from t = 0 to t = T, then solve the 
adjoint system equation, Eq. (A.7), backward starting from t = T with the 
terminal condition, Eq. (A. 8), to t = 0, and at the end of each time step during 
the solution of adjoint system, compute the functional derivative of J^g with 
respect to permeability, Eq. (A. 9), or porosity, Eq. (A.10), at the simulator 
grid cells. Then, one computes the derivative of J^g with respect to the 
spline coefficient W a , Eq. (A. 11), the derivative of Jgf with respect to W a , 
Eqs. (A. 12-19), and the derivative of Jg^ with respect to W a , Eq. (A. 20). 

Because of the large dimensionality of W a o , one seeks to use an 

X, X,y 

algorithm that is as efficient as possible. The essential consideration in the 


choice of a method is the computational time needed to minimize the objective 
function. Most of multivariate minimization methods can be divided into two 
groups: conjugate gradient methods and quasi-Newton methods. The quasi- 
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Newton methods are preferred for moderate sized problems, but the conjugate 
gradient methods become superior to the quasi-Newton methods as the number of 
variables gets large (Scales 17 suggested 250 as a turning point). Although we 
treat 30-201} variables in our examples, the number is larger for field 
applications. The conjugate gradient algorithm of Nazareth 9 was chosen as the 
core minimization method in the present code. 

The remainder of this work is devoted to the numerical evaluation of the 
algorithm on the estimation of permeability and porosity distributions in a 
single-phase, two-dimensional areal reservoir, as modeled by Eq. (2). We want 
to evaluate the algorithm on a well-defined test problem for which the "true" 
property distributions are known a priori. For this reason, we will specify 
the true parameter values, generate our own pressure data by solving the 
reservoir model with these values, and then try to recover the true parameter 
values by using the history matching algorithm. 

Permeability level and distribution is the principal reservoir property 
used to match pressure behavior. Porosity is usually better known than 
permeability, and values from log and core data are often used as initial 
guesses for <j>. (Porosity in the aquifer is generally less well known than in 
the reservoir itself and can be more readily varied than <j> in the reservoir.) 
Aside from aquifer permeability and porosity which are generally not well 
known, reservoir permeability is usually more uncertain than porosity. 

It is difficult to determine the optimal value of the regularization 
parameter even if we know the statistical properties of measurement error of 
the well pressure data. We will choose a set of values of the regularization 
parameter so that they form a geometric sequence and determine the optimal 
regularization parameter from the "quasioptimal" value of the regularization 



13 


parameter 18 which minimizes | |g a 3w a /3$ a | | 2 . 

The data for the cases we will study are given in Table 1 . Although this 
set of data is hypothetical, every effort has been made to have the example 
conform to an actual field simulation. 

An important question concerns starting the algorithm. Convergence 
difficulties are sometimes experienced when the initial guesses of the 
parameters are far from their actual values. To attempt to alleviate this 
problem and to generate an algorithm that is as '’hands-off’' as possible, we 
begin the estimation by determining the unknown parameters as uniform over the 
entire region. Thus, to start, we estimate single values of k and <j> for the 
entire region, called k and <J>, that minimize These values then serve as 

starting points for the full history matching algorithm. The rationale behind 
this strategy is that convergence problems should not be encountered in 
estimating a single parameter. The single value, while not accurate in its 
spatial detail, nevertheless serves as a good starting point for the full 
algorithm. This strategy has been employed in the results to be presented 
shortly. The single variable minimization is carrried out in our code by the 
secant method. 

Table 2 gives the results of this first step for the estimation of <p when 
k is known and the estimation of k when 4> is known. Listed in Table 2 are the 
true values of <{> and k, the initial guesses to start the second method, the 
minimizing parameter value a(oc = <p or k), and the values of Jls» Jst» and J SM 
for various values of the regularization parameter $ at the minimum. The true 
<f> and k surfaces are shown in Figures 2 and 3, respectively. To simulate 
measurement error, uniformly distributed random numbers are added to the 
pressure data generated from our presumed true permeability and porosity 
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distributions. 


The resulting data are shown in Figure *1. 
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ESTIMATION OF SPATIALLY-VARYING PERMEABILITY AND POROSTTY 

We will investigate the effect of the choice of regularization parameter 
(B a ), degree of spline approximation (N xs xNy S ), and the number of observation 
wells in the estimation of <p when k is known and the estimation of k when <J> is 
known. 

We will use six observation wells for values of the regularization 
parameter Bq = 0, 0.0 1, 0.1, 1, and 10 atm 2 for the estimation of <p and S|< = 0, 
0.01, 0.1, 1, and 10 atmVdarcies 2 for the estimation of k and spline grids 
N X sxNys = 5x6, 7x9, and 12x17, where 12x17 is the maximum possible value for 
the pressure grid we are using for both the estimation of <f> and the estimation 
of k. As a special case, we will use 18 observation wells for the estimation 
of k with 8^ = 1 atmVdarcies 2 adn N xs xNy S = 7x9. Finally, in all our 
simulations, we assumed that the root mean square error in pressure 
measurements to be o = 0.3 atm (thus o/p 0 = 0.255 and the corresponding J^g 
value is 18.96 atm 2 ). 

The estimation of <j> started with uniform value of <J> = 0.18*1 and the 
smoothing functional Jg^ was minimized until the change in spline coefficient 
W^, and the gradient of JgM with respect to W^ satisfy the convergence 
criteria given by 

| | W <J>,new _ w <+),old| | < Ei (I4.a) 

||G SM w4) M<e 2 (14. b) 

The same strategy was employed for the estimation of k, where the starting 
value of k = 0.241 darcies (0.243 darcies for 18 observation wells) was used. 
Tables 3 and 4 summarize the history matching results for all the cases 


studied. 
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Effect of Regularization 

We expect that, as the regularization parameter increases, the value of 
the least-squares objective function, Jpg at convergence will increase, but the 
stabilizing functional, Jg-p, will decrease, since a larger value of 
regularization parameter means more regularization on the parameter to be 
estimated at the expense of less exact fitting of the observed well pressure 
data. 

This expectation turned out to be true, with some exceptions for the 
terms JsT^ m+ ^> m = 0 and 1, in the stabilizing functional during the 
estimation of k as shown in Table 4(a). That is, Jg-pd) for 6k = 0.1 
atm 2 / darcies 2 is slightly greater than that for 6k = 0.01 atm 2 /darcies 2 and 
JST^ for Bk = 10 atm 2 / darcies 2 is greater than those for 6k = 0.1 and 1 
atm 2 / darcies 2 . But, the total stabilizing functional, Jgx* and its component, 
Jst^), that represents the third order derivative term and is most important 
among the four terms, Jgx^ m+ ^» m = 0, 1, 2, and 3> decreases strictly without 
exception as 6k increases. 

Table 3(a) shows that Jg?^ terms for true <J> is close to that for 
estimated <{> with 64, = 1 atm 2 for which Jpg and e^JsT are balanced, but Table 
4(a) shows that Jg-pW terms for true k is close to that for estimated k with 
&<p is betweeen 0.1 and 0.01 atm 2 / darcies 2 for which Jpg is 10 - 100 times as 
large as 6k Jst* This means k can be regularized more easily than <]> can. It 
is interesting to compare these numerical indicators of performance with the 
surfaces and profiles in Figures 5 and 6 . The estimated parameter surfaces 
are too bumpy compared to the true surfaces (Figures 2 and 3) when small 
regularization parameters are used and in the non- regularized case. Note the 
"bump" in Figure 5.a or the inflection point at y/y L = 0 .8 in curves 1(6 k =0) and 
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2(6k=0.01 atmVdarcies 2 ) of Figure 6. On the other hand, the parameter 
estimates become too flat for large values of regularization parameters, as 
shown in Figure 5.c and curve 5 in Figure 6, as compared to Figure 2 and curve 
6 in Figure 6. 

To determine the optimal regularization parameter, | |8 a dW a /df$ a | U can be 
approximated by | |(W 2 a -W 1 a )/(£n 6 a> 2 -J ' n 8 a ,1 )| |z where 3a, 1 and 8 a ,2 denote two 

different regularization parameters and Wj 01 and W 2 a , respectively, are the 
corresponding spline coefficients that minimize Jqm and 1 1’| jz denotes Euclidean 
vector norm. Table 5 summarizes the results for the estimation of <)> and the 

estimation of k. The optimal 3<j>= 0.1 to 1 atm 2 for the estimation of <j> and 
the optimal 8k * 0.01 to 0.1 atmVdarcies 2 for the estimation of k which agree 
with the above investigation. 

Table 4 also shows that, in the estimation of k, the Jst^ term, which is 
proportional to the Euclidean norm of k(x,y), in the domain 0 is smaller than 
that calculated from the true k. This behavior can be explained from Eq. (2), 
which shows that the pressure value is governed by the gradient of k with 
respect to the space variables rather than the value of k itself. Thus the 
value of k can be reduced to some extent without changing the values of 
pressure significantly during the estimation of k. 

Effect of Spline Approximation (N xs xNy S ) 

The measurement error for true parameters) is 18.96 atm 2 while 

for B<() = 0 is 20.67 atm 2 for the estimation of <J> and Jls for 8k = 0 is 21.26 
atm 2 for the estimation of k. This can be explained by the fact that the 
spline approximation has the effect of smoothing instead of fitting the noisy 


measured data in detail. 
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It is clear that the measured data can be better fit with more 
parameters, thus we expect that values of Jls should decrease as the dimension 
of spline grid, N xs xNy S , increases. At the same time the estimated parameters 
are expected to be less regular for larger values of N xs xNy S . This 
expectation turns out to be true in our examples in the estimation of both <p 
and k. The value of JgT^ for the true <f> is closer to that for the estimated 
<t> when N xs xNy S = 7x9, and the values of JgT^ for the true k lies between 
those of the estimated k for N xs xNy S = 7x9 and 12x17. However, if we want the 
values of JgT^ n °t to be greater than that of true k, we conclude that 
N xs xNys = 7x9 is the best value for the 10x15 pressure grid in our example. 

Figures 7 and 8 show the effect of spline approximation on the estimation 
of <p and k, respectively. We can observe ill- conditioning in the estimation of 
k from Figure 8c (N xs xNy S = 12x17) which can be explained since the ratio of 
spline grid to pressure grid (h) is only 1.09, so that we have 204 unknown 
spline coefficients to be estimated as compared with 210 measured pressure 
data. 

Effect of Number of Observation Wells 

The regularization effect is relatively less important for more 
observation wells and thus the values of Jg? and the terms in Jg-p, Jg T ( ra+ ' 1 ), m 
= 0, 1, 2, and 3 are larger as shown in Table 4(c). As one can see in Figure 
9, the estimated k for 18 observation wells is closer to the true k than that 
for six observation wells. Note that J^g for 18 observation wells case is 
closer to the true value than that for six observation wells with any set of k 
and N xs xN ys is. 



19 


CONVERGENCE OF THE ALGORITHM 

When we seek the minimum of Jqm with respect to W a (ot = <J> or k) with the 
conjugate gradient algorithm, we need up to N xs xNy S conjugate directions to 
find the approximation of the inverse of the Hessian matrix. However, the 
algorithm employed here uses 5 to 10 conjugate directions to find the 
approximation of the inverse of the Hessian matrix. On the whole, the 
algorithm requires 1 0 to 20 different inverse Hessian matrix evaluations which 
correspond to 0.5 to 1 hour of CPU time on a VAX1 1-780 for a single history 
match. 

The algorithm determines the minimum of JsM(W a + s6W a ) along s, the step 
size, by trial and error. In some cases, if s is too large, some elements of 
W a + sSW a have negative values that are physically impossible. Thus, in the 
implementation of the algorithm a limit on the size of s was employed. 
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CONCLUSIONS 

In this study we have developed and tested an automatic history matching 
algorithm for estimating spatially-varying porosity and permeability in a 
single phase areal reservoir. The algorithm is based on spline approximations 
of the parameters and a regularization formulation. In the regularization 
approach to parameter estimation by introducing the stabilizing functional as a 
measure of "non- smoothness,” one can control the properties of the para m eter 
estimates as well as the history match. 

We have presented results of a detailed numerical evaluation of the 
performance of the history matching method. In this example, we found that 
the permeability distribution is estimated somewhat better than the porosity 
distribution at comparable levels of spline approximation and degree of 
regularization. It was also found that increasing the value of the 
regularization parameter leads to estimated property distributions that are 
smoother than those obtained for smaller values of the regularization 
parameter. Some exceptions to this behavior were found at small values of the 
regularization parameter that can be attributed to inherent numerical ill- 
conditioning in estimation problems of this size. There appears to be an 
optimal level of spline approximation in the case of the example studied. The 
optimum was a 7x9 grid, for which the ratio of the size of the spline grid to 
that of the pressure grid is 2.5. This optimal value of approximation appears 
to represent a trade-off between a low dimensional spline grid that has as few 
unknown parameters as possible and a high dimensional spline grid that is more 
able to represent the details of property distributions but introduces more 
unknowns and therefore inherently more ill- conditioning in the optimization 
step. 
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Based upon the optimal spline approximation, the optimal regularization 
parameters are - 0.1 to 1 atm 2 for the estimation of <j> and Bk = 0.01 to 0.1 
atmVdarcies 2 for the estimation of k which are determined from the 
"quasioptimal" condition of regularization and give the same magnitude of the 
values of measure of "non-smoothness" compared to the true profiles. 

Finally, we can suggest a history matching strategy: 

1. Choose simulator and spline grid systems. The number at the spline 
coefficients need not be as large as the number of simulator grid 
cells. 

2. Find uniform initial guess of parameter to be estimated which minimizes 
Jls and calculate Jls/^ST at convergence. 

3. Choose the regularization parameter value approximately the same as 
J LS/JST above and find a set of spline coefficients that minimizes JgM* 
Step (3) can be repeated to evaluate the result for the different 
regularization parameter values around the J^s/^ST value determined in 
Step (2), so that we can find the "optimum" value of regularization 
parameter discussed in the previous section. 
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APPENDIX A. GRADIENT OF THE OBJECTIVE FUNCTION WITH RESPECT TO THE UNKNOWN 
PARAMETERS 


Let us begin by supposing that we want to minimize the least-squares 
objective function, J^^, given by Eq. (5), with respect to permeability k 
and porosity <p subject to Eqs. (2)-(4). By adjoining the model of Eq. (2) 
to J^ by means of an adjoint function i^(x,y,t). 



+ Mx.y.t) [ -c<t»-j^p(x,y,t) + v«(£ V p(x,y,t)) 


W=1 


dtdxdy 


(A.l) 


then minimizing J^<. leads to the following equations governing ^(x,y,t), 
c<P iMx.y.t) = - ^(x,y,t)) 


{l p M- 

n=l k=l 

p(x,y,t) ) 5(x-x k )5(y-y k )6(t-t n ) 
in ft x [0,T] 

} 

(A. 2) 

|| = ° 
dn 

on 8ft x [0,T] 

(A. 3) 

V (x,y,T) = 0 

in ft 

(A. 4) 


where n has the direction outward normal to the boundary, and the functional 
derivatives of J. ~ with respect to k and cf> at (x,y) 6 ft are 
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(x,y) = - 



Vij;(x,y,t)-Vp(x,y,t) dt 
in ft 


(A. 5) 


G^ s (x,y) = - c j 4>(x,y,t) p(x,y,t) dt 
o 

/•T 

■ ciJ;(x,y J 0)p o (x,y) + c / (» i|>(x,y,t)J p(x,y,t) j dt 

° in « (A. 6) 


The first order necessary condition for a local minimum of J LS is that p 

and i p satisfy Eqs. (2)-(4) and (A.2)-(A.4), respectively, and that 

G^$(x,y) = 0 (for the estimation of k) or G^(x,y) = 0 (for the estimation of <j>) 

k 

for all (x,y) € ft. The gradients G^<-(x,y) and G^(x,y) are used in the so- 

called optimal control algorithms for history matching. As noted, since these 
gradients can be calculated directly without requiring the sensitivity coefficients, 
3p/3k and 3p/34>, these optimal control algorithms are computationally attractive 
for history matching. 

The adjoint equations (A.2)-(A.4) can be written in a finite difference 
form corresponding to Eq. (12) as 


- ♦?> "E - *1 J 

j€J i 

N obs 

-2^ (p? bs - p n )S. • 
L-+ H k,n 1 ,1. 

k=l 

for i ( and n = 1,2,..., 


N.+l 

V =0 1 € \ 


(A. 7) 


(A. 8) 



24 


where J. = { i -N x , i-1, i+1, i+N x > nlN^ and the derivatives of with 
respect to and (J> i are 


n=l j€Jj 


(A. 9) 


. = 

b LS ' 


N t for 1 £ 

n, n _n-l- 




n=l 


N 


= Qc { ^ 

n=l 


(A. 10) 


for i € W. 


In our algorithm, k and <p are represented by the bi-cubic spline approxi- 
mation Eq. (8), and the actual unknown parameters are the coefficients, 


W 


£ , £ 

x’ y 


and W 


<t> 

vy 


Thus, we need to obtain the derivatives of the overall 


k 6 

objective functional with respect to W. . and WT . . These gradients 

x’ y x’ y 


are then the values used directly in the conjugate gradient minimization 
method. 

Let the N x N matrix B have elements B . • = b (£ ,(i - -ijAx) 

XS X X X 5 )C j 1 XX X L 

XX i 

and the N x N matrix B have elements B . . = b (£ ,(i - -i)Ay). Then 

ys y y y»~y**y y y y ^ 

the derivative of with respect to the elements of W a (a = k or <J>) is 


..a 


LS 


y 

,£ v »£ u " 2_i Lu B x,£ x , i x G LS,i B y,£ w ,i 


x' y 


v 1 y 1 


y y 


£ x = 1,..., N xs and £ y = 1,..., N y$ (A. 11) 

where i = i + N (i -1). Thus Eq. (A. 11) relates the gradient of the least- 
x x y 

squares objective function with respect to the spline coefficients to that 
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with respect to the individual reservoir parameters at each grid point of 
the simulator. 

Eq. (A. 11) expresses the gradient of the least-squares portion, 
of the overall objective function, We need to obtain the gradient of 

Let us consider the second component of namely the stabilizing 
functional, J ST , given by Eq. (6) which is 

J = r + r j(2) . r j(3) , r i(4) /. \ 

J ST C 1 J ST ? 2 d ST c 3 d ST + c 4 d ST (A.12.1) 


N N 



(A. 12. 2) 


(A. 12.3) 


(A. 12. 4) 


M N 



(A. 12. 5) 


where 4 ^ >0, ^ _> 0, ^ >. 0, 4 ^ > 0 and E, = x/Ax and n = y/Ay. Using E, and n 
rather than x and y, b x U x ,x) and b U ,y) in Eqs. (9) and (10) and their 
derivatives are 

b x (V x ) = X * 4 ( 4 ^x + ? /h x ) (A. 13.1) 

b y U y ,y) = X* 4 (4-£ y + n/h y ) (A. 13.2) 

^5 b x ( V x) ■ X* 4(m >(4-i x + c/h x ), m = 1,2,3 

X 


(A. 13.3) 
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Tli W y) = IT X‘ 4(m) (4-t v * n/h v ), m = 1,2,3 


dn 


where x* 4 ^(*) denotes m-th derivatives of X* 4 (*)> h = Ax /Ax, 

X s 

Ay /Ay. Combining Eqs. (8), (A. 12), and (A. 13) 


where 


u _ o fl 0x -Oy 

ST,£ >£ ,m ,m 1 £ ,m £ , 


x y 7 x’ y 


x' x y' m y 


+ c (a 1x A 0y + A^ x A^ y ) 

2 l V m x VV V”x y m y / 

+ £ (a 2x A°y + 2A^ X 

3 \ £ x’ m x £ y’ m y £ x’ m x £ y’ m y 


+ A? X A 3y 
£ x’ m x £ y’ m y 


+ C (a 2x A 0 ^ + 3A 2x a^ 

4 \ £ x ,m x £ y’ m y £ x’ m x £ y’ 

+ 3A* X A 2 ^ + A^ x A 2 ^ \ 

£ x’ m x £ y’ m y £ x’ m x £ y’ m y) 


y ,m y 


,w 


a 


N N 
xs ys 


G ST,l v ,* * 2_! 2-. H ST,1 ,1 ,m„,- " 


,a 


y m t =1 m =1 
x y 

N N 
xs vs 


x’ y’ x’Hy ^Vy 


ST 


- V 4 Vlr w Ii a 
Lu ST,£ ,£ w £ , 


£ =1 £ =1 
X y 


"xs' 3 


X' y x’ y 




(A. 13.4) 
and h^ = 


(A. 14) 
(A. 15) 
(A. 16) 


m = 0,1, 2, 3 


(A. 17) 
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A 


My 

V £ y 









m = 0,1, 2, 3 ( A . 18) 


or in the matrix form. 


A 0X = 



A 


lx 



20 

129 

60 

1 




129 

1208 

1062 

120 

1 



60 

1062 

2396 

1191 

120 

1 


1 

120 

1191 

2416 

1191 

120 

1 


1 120 

1191 

2416 

1191 

120 

1 

1 

120 

1191 

2396 

1062 

60 


1 

120 

1062 

1208 

129 



1 

60 

129 

20 


6 7 -12 -1 


7 

40 

-22 

-24 

-1 





-12 

-22 

74 

-15 

-24 

-1 




-1 

-24 

-15 

80 

-15 

-24 

-1 





-1 

1 *3" 

1 CM 

1 1 

-15 

80 

-15 

-24 

-1 




-1 

-24 

-15 

74 

-22 

-12 





-1 

-24 

-22 

40 

7 






-1 

-12 

7 

6 
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(A. 19) 

and A 0- ^, A^, A^, and A 3 ^ have the analogous expressions to Eq. (A. 19). 

Finally, we obtain the gradient of overall objective function J_,. 

Src 

with respect to W a as 

r W a _ r w a W a 

b SM,£ x £ y " G LS,£ x £ y 3 a G ST,£ x £ y 


(A. 20) 
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Weighting Factor c m+1 > m = 0, 1, 2, and 3 in Eq. (A. 12.1) 


In the stabilizing functional in Eq. (A. 12.1), we have four undetermined 

constants, namely the weighting factors, £ m = 0, 1, 2, and 3, which are 

arbitrary except for the conditions that > 0, > 0, ^ ^ 0, and > U. 

We can set these constants in a systematic way by using the fact that we want 

the four terms in Eq. (A. 12.1), each of which is a weighting factor 

multiplied by J^ + ^, described in Eqs. (A. 12.2)-(A.12.5) , to be of about 

equal magnitude. We assume that h , the ratio of the size of the spline grid 

to that of the simulator grid along x-direction is not much different from hy, 

the ratio along the y-direction, and we let h = (h h ) 2 . Then, the ratio of 

x y 

terms J^y , m = 0, 1, 2, and 3 in Eq. (A. 12) is approximately 



and this suggests values for the weighting factors, 5 y, m = 0, 1, 2, and 3 
2 4 6 

as = 1, ^2 = h , £3 = h , = h , where is independent of h. 




30 


A (m)x 

H £ m 
X X 

My 

a £ m 

y y 

b x ( V x) 

b y (^y.y) 

B.. 


NOTATION 

matrix defined by Eq. (A. 19), m = 0,1,2, and 3 

matrix defined by Eq. (A. 19), m = 0,1,2, and 3 

quantity defined by Eq. (A. 17) 

quantity defined by Eq. (A. 18) 

cubic spline function defined by Eq. (9) 

cubic spline function defined by Eq. (10) 

N xN matrix of spline function values 

A J X 

N ys xNy matrix of spline function values 
compressibility, atm” 1 [Pa *] 


G^ s (x,y), G^ s (x,y) functional derivative of J LS with respect to k 
and <p at (x,y) 

derivative of with respect to the values of k and <f> 
at the grid point of the simulator, i = 1,2,..., N 

w k ^4> 

\s £ £ G LS £ 9 derivative of J LS res P ect to the values of 


pk 

G LS,i’ G LS,i 


W £ i and W * i 
x * y x * y 


,W 


,w 


<t> 


SM £ £ U SM £ £ derivative of with respect to the values of 

x y x y 


w i i and i 
x’ y x 9 y 

G ST £ £ » G 5 j £ £ derivative of with respect to the values of 

x y x y 


w t . and W* 
x ,X/ y x’y 
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H 3 (£!) 


ST,£ x ,5, y ,m x ,m y 

h 

h . 


’’Vy 


'LS 


J SM 

ST 

(m+1) 

ST 


J . 
1 


k i 


Sobolev space of order 3 on the domain Q, 
quantity defined by Eq. (A. 14) 

< h xV* 

AX g /AX 

Ay s /Ay 

indices for simulator grid, i = i + N (i -1) 
least squares objective function, Eq. (5) 
smoothing functional, Eq. (7) 
stabilizing functional, Eq. (6) 

terms in stabilizing functional defined by Eq. (A. 12), 
m = 0,1,2, and 3 

a set of integers which indicate the neighborhood of 
i-th grid block 

2 

permeability, darcies [m ] 

2 

permeability values at thei-th simulator grid, darcies [m ], 


i = 1 


9 


N 


V £ y 


obs 


’w 


My 

N xs’ N ys 

V 


indices for bi-cubic spline approximation grid, £ = 1,..., 

X 

N xs’ £ y = 1 N ys 

N N , total number of simulator grid blocks 
x y 

number of observation locations 

number of wells 

number of observation times 

number of simulator grid blocks along x- and y-directions 
number of nodes along x- and y-direction of spline grid 
set of integers from 1 to N 
set of integers from 1 to N x and N y 
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pressure, atm [Pa] 
observed pressure, atm [Pa] 
initial pressure, atm [Pa] 

volumetric flow rate per thickness of reservoir of well w, ft^/day 
[m /s] 

AxAy/At 

Ay/ Ax or Ax/Ay 
time, days [s] 

time periods over which observations are available, days [s] 

N xs xN ys matnx of spline coefficient of k and <f> 

spatial variable, miles [m] 

extent of domain in x-direction, miles [m] 

spatial variable, miles [m] 

extent of domain in y-di recti on, miles [m] 


Greek Letters 

a unknown parameter to be estimated (a = k or <}>) 

&l< regularization parameter for the estimation of k, 

atm^/darcies^ [Pa^/m^] 

^ regularization parameter for the estimation of <j>, atm^ [^ a ] 

At time interval of observation, days [s] 

AX, Ay simulator grid spacings, mile M 

AX s ,Ay $ spline grid spacings, mile [m] 

6(*) Dirac delta function 



? m+l 


n 

u 


Kronecker delta 

weighting factor of Sobolev norm 

y/ Ay 

viscosity, centipoise [P a *s] 


INI o , m = 0, 1 , 2, and 3 
rfa) 
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£ X/AX 

a standard deviation of measurement error, atm 

4> porosity 

4 

X* (•) cubic B-spline function, Eq. (11) 

i p adjoint variable 

ft spatial domain of reservoir 

3ft boundary of reservoir 

SI Metric Conversion Factors 


* 


atm 

xl. 01325 

E+05 

= 

Pa 

2 2 
atm /darcies 

xl .054 

E+34 


Pa 2 /m 

cent i poise 

xl* 

E-03 

= 

Pa.s 

darcies 

X9.869 

E-13 

= 

m 2 

days 

xl . 157 

E-05 

= 

s 

ft 2 /day 

xl .075 

E-0.6 

= 

2. 
m /s 

miles 

xl. 609344* 

E+03 

= 

m 


Conversion factor is exact. 
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Table 1 . Specifications of Reservoir for History Matching Example 



Oilfield units 

SI units 

Dimension of reservoir 

12.4x1 8.6(mi) 

20, 000x30, 000(m) 

Compressibility of system 

1.2x10“ 5 (atm~ 1 ) 

1 .2x10~ 10 (Pa“ I ) 

Viscosity of fluid 

2.0 (cp) 

2.0x1 0' 3 (Pa.s) 

Number of production wells 

1 

1 

Production rate 

500 (ft 2 / day) 

5.376x1 0‘ 4 (m7s) 

Initial pressure 

150(atm) 

1.52x10 7 (Pa) 

Pressure grid 

10x15 

10x15 

Number of observation wells 

6 

6 

Number of well pressure data per each well 35 

35 

Time interval of well pressure data 

10 (days) 

8.64x1 0 6 (s) 

Total number of data points 

210 

210 
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Table 2. True Values of Porosity arid Permeability, the Uniform 

Values of <f> and k that Minimize and the Corresponding 
Starting Values of J L <., J<.y, and 


Parameter to be 
estimated 

4> 

k' 1 ’ 

True Value 

0.2-0.05 sin(irx/x^) 

0. 3-0.1 sin (ttx/x^) 


•sin (2 uy/y L ) 

•sin (27iy/y L ) 

Initial guess 

0.25 and 0.15 

0.25 and 0.35 

a that minimizes 
J LS 

0.184 

0.241 

J LS 

44.6 

41.0 

°ST 

5.1 

8.7 


o < 2) 

44.6 

. 0 1 2 (3) 

41 .0 

0.01 

44.7 

0.01 

41.1 

0.1 

45.1 

0.1 

41 .9 

1 

49.7 

1 

49.7 

10 

95.6 

10 

128.0 


(1) units are darcies 

(2) units are atm 

2 2 

(3) units are atm /darcies 
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Table 3. Estimation of 4> 

(a) Final values of performance indices for o = 0.3 atm and 
N xs xN ys = 7x9 as a fL,nctl ’ on of B^ (atnr) 



J SM 

a 

J LS J ST 

J (1) 

J ST 

jffxlO 

2 j(3>Kl0 3 

J^xlO 3 

0 

20.67 

20.67 15.55 

6.47 

5.30 

30.24 

193.63 

0.01 

20.86 

20.76 9.47 

6.47 

2.67 

10.96 

61 .64 

0.1 

21 .60 

20.85 7.48 

6.44 

1.76 

5.11 

18.75 

1 

27.61 

21.16 6.45 

6.28 

0.75 

1.47 

1.79 

10 

81 .69 

24.09 5.76 

5.74 

0.20 

0.12 

0.09 

true 


7.09 

6.17 

2.63 

7.48 

1 .88 

V> 

C 2 = h 2 = 

» 4 3 =h 4 , C 4 =h 6 , with 

h = 2.5 




(b) Final values of performance indices for a = 0.3 

(3, = 1 atm^ as a function of N xN 
4> xs ys 

atm and 


N xN 
xs ys 

h 

a 

J SM J LS J ST 

j(l) 

J ST 



J^>xl0 3 

5x6 

5 

27.86 21.38 6.47 

6.29 

0.50 

0.04 

0.002 

7x9 

2.5 

27.61 21.16 6.45 

6.28 

0.75 

1.47 

1.79 

12x17 

1 .09 

27.17 20.78 6.40 

6.31 

3.11 

18.38 

14.78 


a q=l, ? 2 =h 2 , C 3 =h 4 , ? 4 =h 6 
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Table 4. Estimation of k 


(a) Final values of performance indices for a = O^atm and 

N xN = 7x9 as a function of &. (atm /darcies ) D 
xs ys k 


6 k 

°SM 

J LS 'Si' 

j(D 

J ST 

0<2>xl0 : 

2 J s?’ x1 ° : 

3 o^xio 3 

0 

21.26 

21.26 20.46 

11.68 

7.81 

36.75 

28.05 

0.01 

21.59 

21.44 14.86 

11.36 

5.56 

18.14 

9.99 

0.1 

22.82 

21.49 13.36 

11.38 

4.79 

12.33 

4.84 

1 

33.42 

22.18 11.24 

10.54 

3.83 

5.36 

1.04 

10 

114.7 

35.70 7.90 

7.39 

4.98 

2.86 

0.33 

true 


17.54 

13.88 

10.50 

29.97 

7.50 

a i 

c r lf 

r -h2 

h , 

C 3 =h^, c^h 6 with h = 

= 2.5 




(b) Final values of performance indices for a = 0.3 atm and 
6^ = 1 atm^/darcies^ as a function of N xs x Ny S 


"xs xN ys 

h 

J SM J LS J ST 9 

j(D 

J ST 



j' 4) xl0 3 

5x6 

5 

35.70 23.48 12.22 

11.30 

1.30 

0.39 

0.02 

7x9 

2.5 

33.42 22.18 11.24 

10.54 

3.83 

5.36 

1.04 

12x17 

1.09 

31.59 21.68 9.91 

9.53 

12.43 

81.60 

68.91 

a V 1 - 

t 2 =h2 > 

{ 3 =h 4 , ? 4 =h 6 . 





(c) Final values of performance indices for o = 0.3 atm, 3^ 

1 atm^/darciesS and N xs xNyS = 7x9 with different number 
observation wells 

of 

# of 
wells 

J SM 

J LS V ‘ 

,(1) j 

] ST J 

l^xlO 2 J 

s ? )xl ° 3 J 

< 4 W 

ST xlu 

6 

33.42 


10.54 

3.83 

5.36 

1.04 

18 

76.96 


11.96 

5.00 

11.28 

3.74 

V- 

^ 2 =h2 » 5 

3 =h^, ^ 4 =h® with h = 

2.5. 





b 7 7 . " . 

To convert from atm /darcies to SI units see Notation. 
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Table 5. Values of || B C || 

a 

(a) Estimation of <f> 




II 6* 

0.01, 

0.1 

0.167 

0.1, 

1 

0.148 

1, 

10 

0.169 


(b) Estimation of k 



*k 

II e k dw k /dB k || / 

0.01, 

0.1 

0.211 

0.1 , 

1 

0.229 

1 , 

10 

0.456 


a Units are darcies- To convert to SI units see Notation. 
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Figure 1. 

Figure 2. 
Figure 3. 
Figure 4. 


Figure 5. 
Figure 6. 


Figure 7. 
Figure 8, 


FIGURE CAPTIONS 

Pressure and spline grid system, pressure grid: 10x15; 
spline grid: 5x6; 7x9; 12x17 

(0 = Observation well; P = production well) 

True <f>(x,y) surface 

(j>(x,y) = 0.2-0.05 sin (trx/x^sin (ETry/y^) 

True k(x,y) surface 

k(x,y) = 0. 3-0.1 sin (iTX/x^)sin(2Try/y L ) darcies 
Simulated pressure data vs. time for a = 0.3 atm 

1. (3.1,3.1) a 

2. (9. 3,3.1) 

3. (3. 1,9. 3) 

4. (9. 3,9. 3) 

5. (3.1,15.5) 

6. (9.3,15.5) 

a units are miles 


Estimated 4> surface for o = 0.3 atm, N xs xN ys 7x9, and 
from top down: 2 

B = 0 atnr 

4> 

0.01 

0.1 

1 

10 


Cross-sectional plot of k(* L /2,y) vs. y for o = 0.3 atm 

"xs xN ys = 7x9 a " d 

2.2 

1. B k = 0 atm /darcies 

2 . 0.01 

3 . 0.1 

4. 1 

5 . 10 

6. true values 

7. initial guess 

Cross-sectional plot of <J>(x L /2,y) vs. y for a = 0.3 atm 

S = 1 atm , and 

4> 


1. N xs xN ys - 5x6 

2. 7x9 

3. 12x17 

4. true values 

5. initial guess 

Estimated k surface for o = 0.3 atm, B k = 1 atm^/darcies and 
from top down: 

N xs xN ys = 5x6 
7x9 
12x17 
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Figure 9. Estimated and true k surfaces for a = 0.3 atm, 

6k = 1 atrrr/darcies , N xN = 7x9, and from top down 
* xs ys 

6 observation wells 
18 observation wells 
true k 
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